In [1]:
%matplotlib notebook
from astropy.cosmology import LambdaCDM
import numpy as np
import matplotlib.pyplot as plt
from astropy import constants as const
import astropy.units as u
from scipy.integrate import quad

In [2]:
cosmo = LambdaCDM(H0=70, Om0=0.3, Ode0=0.7, Tcmb0=2.725)

In [3]:
z = np.arange(0, 1.5, 0.01)
z=0.4

In [4]:
phi_star = 3.6 * cosmo.efunc(z)**2
alpha = -1.05 * (1 + z)**(-2/3)
fr = 0.8*(1 + z)**(-1/2)

In [5]:
alpha


Out[5]:
-0.8390167065610477

Need to compute the cluster volume...

$M_{vir} = 4/3 \pi r^3_{vir} \rho_c(r<r_{vir}) = 4/3 \pi r^3_{vir} \Delta_c \rho_c$

if we let $\Delta_c = 200$ then

$M_{200} = 4/3 \pi r^3_{200} 200 \rho_c$ with $\rho_c = \frac{3H(z)^2}{8\pi G}$

or just $M_{200} = V_{200}200\rho_c$

Don't forget that $H(z) = H_0E(z)$


In [6]:
def rho_crit(z, cosmo):
    # convert G into better units:
    G = const.G.to(u.km**2 * u.Mpc/(u.M_sun * u.s**2))
    return 3 / (8 * np.pi * G) * cosmo.H0**2 * cosmo.efunc(z)**2 # Mpc^3

# So now we are going to calculate the volumes as a function of z
M200 = 1e15 * u.solMass
V200 = M200/ (200 * rho_crit(z, cosmo))

In [7]:
V200


Out[7]:
$24.135203 \; \mathrm{Mpc^{3}}$

The Schechter Function:

For Luminosity:

$\Phi(L) = \phi^\star \frac{L}{L_\star}^\alpha e^{-\frac{L}{L_\star}}$

For Magnitudes:

$\Phi(M) = \phi^\star\frac{2}{5}log(10) (10^{\frac{2}{5}(M_\star - M)})^{\alpha+1} e^{-10^{\frac{2}{5}(M_\star - M)}}$


In [8]:
def schechterL(luminosity, phiStar, alpha, LStar): 
    """Schechter luminosity function.""" 
    LOverLStar = (luminosity/LStar) 
    return (phiStar/LStar) * LOverLStar**alpha * np.exp(- LOverLStar) 

def schechterM(magnitude, phiStar, alpha, MStar): 
    """Schechter luminosity function by magnitudes.""" 
    MStarMinM = 0.4 * (MStar - magnitude) 
    try:
        return (0.4 * np.log(10) * phiStar * 10.0**(MStarMinM * (alpha + 1.)) * np.exp(-10.**MStarMinM))
    except OverflowError:
        return 0.0

In [9]:
M200 = 10 # 10^14 Msun
Mpiv = 6 # 10^14
zpiv = 0.6

alpha = -0.96 * (M200 / Mpiv)**0.01 * ((1 + z)/ (1 + zpiv))**-0.94
Phi = 1.68 * (M200 / Mpiv)**0.09 * ((1 + z)/ (1 + zpiv))**0.09 * cosmo.efunc(z)**2
fr = 0.62 * (M200 / Mpiv)**0.08 * ((1 + z)/ (1 + zpiv))** -0.80

In [10]:
M_star = -20.44 + 5 * np.log10(0.7) # abs mag.
M_star = -21.946521706412092 # abs mag @ z= 0.4
M_star_sub = M_star - 2.5 * np.log10(0.4)
M_star_big = M_star - 2.5 * np.log10(100)
#M_star = 20.071802735723015 # app mag @ z= 0.4
y, err = quad(schechterM, -np.inf, -21.61828907152557, args=(phi_star, alpha, M_star))
#y, err = quad(schechterM, -np.inf, M_star_sub, args=(phi_star, alpha, M_star))

In [11]:
(y * V200.value + 1) * fr


Out[11]:
33.090792933096466

In [12]:
Phi


Out[12]:
2.6476407471138423

In [13]:
phi_star


Out[13]:
5.484097445491274

In [14]:
1/5


Out[14]:
0.2

Here is the stuff from the c++ code.


In [22]:
def dAintegrand2(y, alpha, mstar):
    #val = y**alpha * np.exp(-y)
    # using mags instead of luminosity
    val = 10.0**(-2 / 5 * (y - mstar) * (alpha + 1.)) * np.exp(-10.** (-2 / 5) * (y - mstar))
    return val
    
def calc_hon(z, mlimit, mstar):
    M200 = 10 # 10^14 Msun
    Mpiv = 6 # 10^14
    zpiv = 0.6
    
    alpha = -0.96 * (M200 / Mpiv)**0.01 * ((1 + z)/ (1 + zpiv))**-0.94
    Phi = 1.68 * (M200 / Mpiv)**0.09 * ((1 + z)/ (1 + zpiv))**0.09 * cosmo.efunc(z)**2
    fr = 0.62 * (M200 / Mpiv)**0.08 * ((1 + z)/ (1 + zpiv))** -0.80
    
    
    R200 = M200 * 4301.8644383/ (cosmo.H0.value**2 * cosmo.efunc(z)**2)
    R200 = R200**(1 / 3)
    
    # Calculate factors
    Vol = 4 / 3 * np.pi * R200**3 # This gives basically the same number as my calculation above.
    # Using Mags Instead of Luminosity
    NN = -0.4 * np.log(10) * Vol * Phi * quad(dAintegrand2, -np.inf, mlimit, args=(alpha, mstar))[0]
    print(NN)
    return (NN + 1) * fr

In [23]:
calc_hon(z, -21, M_star)


-inf
/home/boada/.local/lib/python3.6/site-packages/ipykernel_launcher.py:4: RuntimeWarning: overflow encountered in exp
  after removing the cwd from sys.path.
Out[23]:
-inf

In [17]:
type(np.log(10))


Out[17]:
numpy.float64

In [18]:
cosmo.efunc(z)**2


Out[18]:
1.523360401525354

In [19]:
from scipy import special

In [20]:
phi_star * special.gammainc(alpha + 1, 3) * V200 * fr


Out[20]:
${\rm NaN} \; \mathrm{Mpc^{3}}$

In [21]:
alpha


Out[21]:
-1.093961783729938

In [ ]: